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PATENT 

INSTITUT FRAN^AIS DU PETROLE 

METHOD FOR MODEtUNG FLUID FLOWS 
IN A FRACTU RED MULTILAYER POROUS MEDIUM 
AND CORRELATIVE INTERACTIONS IN A PRODUCTION WELL 

Inventor : .Sylvain SARDA 
ABSTRACT 

Me^^e dA method is disclosed for modelling fluid flows in a fractured multilayer 
porous medium by taldng aooouiit o f accounting for (he real geometry of the fracture 
network and the local exchanges between the porous matrix and the fiactures at each 
node of the network, thus allowing simulation of the interactions between the pressure 
and flow rate variations in a well running across said^ meditim. The method 
essentially comprises disoretiggtion o f discretizing the fractured medium by means of a 
mesh pattern, with fracture meshes ocntr e d centered on nodes at the various intersections 
of the fracturesi with each node being associated . with a matrix volume, and 
determination of the flows between each fracture mesh and the associated matrix 
volume in a pseudosteady state. The method can be applied in hydrocarbon production 
well testing for ojcomplo . 
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BACKGROUND OF THE INVENTION 
KfF . i , n art TITF. IN\fBNTION FIeid of the Invention 

The present invention relates to a method for simulating well tests on a discrete 
fractured reservoir model. 

5 ¥ho method can bo implem e nt e d for exampl e in ih & fiold of oil production by 
roaorvoir onginooro in order to obtain reliabl e flow prodiotiono. 

A - well tost iB to bo aimulat e d in a pomioablo porous modium (res e rvoir) cro o cod 
through by a network of fraotuico - much more oonduoting than tho poroua - matrix. 
Fraotared rooervoirg oonotituto an oxtr e mo t>Tio of hoterog e nooup looorvoiro oompriisin g 
10 two v e ry different media, a matrix medium containing tho moat port of tho oil in plaoo 
and having a low permeability, and a fractur e d modium repr e s e nting loflo than 1 V o of 
th e oil in plac e and highly conducting. Tho fractur e d m e dium itaolf can bo complex, 
with d i fferent o e ri ois of f ru e turc s charactorizod by - their r e spoctivo donaity, length^ 
orientation, inclination and op e ning. 

15 During woU tooting, the flow rate conditions impoaod on tho well l e ad th e oil 

co ntain e d in tho rosorvoir to flow towards tho well. It i s a singl e phase (the otHy mobil e 
phacp is tho oil) and comproGGiblo flow (the rock of tho r oc o f veir ond the oil fluid ar e 
oompres s ibl e )- 

FoT any olomontoty volum e of tho rooor^'oiTj the pr e ssuro of tfao oil contained in this 
20 volume ia govern e d by tho following eqxiation, if gxavity la not tak e n into account : 



PAGE 2(146 ' RCVD AT nil4 3:51:34 PM (Eastern Standard TUne] ' SVR:USPT0{FXRF-2lj ' DNI5:7469166' CSID:703 312 6666 ' DURATION |inni-ss);10-18 



03/03/04 WED 16:57 FAX 703 312 8666 A T S K ©027 



^ : rcpreaontQ the pore volume ; 
Co:, tho total comproooibiUty (fluid + rock), 
K, th e ppimcability of the rock, 
5 \iy tho viooosity of the fluid, 

Q, tfao inooming flow, whioh is aoro Gveryv ^ ^horo oxoopt in plaooo whoro tho woll 
oommimiootoa with tho roaon^oir^ and 
P, tho unknown proQOur e . 

la ord e r to oimulate ft woll teat, whatovor - thc m e dium, thio e qimtion hno to h e 
10 golvod in opaoo and in timo. Digorotiaation of tho ro6or>^oir (moah pattom) iq thoroforc 
pcrform o d and o olution of tho probl e m conoiato in finding tho prossuros of tho m e sh e ti 
with tim e , its e lf diGorotiaod in a oottain numb e r of tim e int e rN^ale r 

In the ooo e of a fraotur o d m o dium, if all th e comploxity of tho ftacturo network is 
to bo takon into account, it io ncoossaiy to have a mcoh pattom that ropregents th ia 
15 notv ^ ^ork aooumtoly. Furthormoro^ in auoh a medium, llio pormoabilitiQii K oro gen e rally 
muoh high e r in flio fraotUTCO than in tho - matrix> and tho flowo ore oomoquontly foot in 
tho fraotutoa and glow in tho matrix. 

BArKf^nnuNP OF THE INVENT lO N DescrlDtion of the Prior Art 

Single-phase flow computing tools are currently use d: hQ> ^ 'Ovor . However , 
20 the vthese tools are not applied to the real geologic medium in ftH-its view of the 
complexity ha tthereof and instead are applied to a homogenized representation, 
according to the reservoir model called a "double-medium moder described for 
example by Wanren and Root in "The Behavior of Naturally Fractured Reservoirs", SPE 



PAGE 27/46 * RCVD AT W004 3:51 :34 PM pastern Standard Time] * SVR:USPT0€FXRF-2/1 * DNIS:F469165 * CSID:703 312 6666 * DURATION (mm-ss):10.18 



03/03/04 WED 16:57 FAX 703 312 6666 A T S K 21028 



Journal, September 1963. Any elementary volume of the fractured reservoir is thus 
n\odelled in the form of a set of identical parallelepipedic blocks limited by an 
Orthogonal system of continuous uniform fractures oriented in the direction of one of 
the three main directions of flow (model referred to as "sugar box" model). The fluid 
5 flow on the reservoir scale occurs through the fractured medium only, and fluid 
exchanges take place locally between the fractures. and-feeThe matrix blocks are not 
applied to the complex medium. This representation, which does not reproduce the 
complexity of the fractuj« network in a reservoir, is however effective but at the level of 
a reservoir mesh whose typical dimensions are 100 m x 100 m. 

10 It is however preferable to carry out the flow calculations on the "real" geologic 
model instead of an equivalent homogeneous model, which affords the double 
advantage of allowing : 

to validate v alidation of the geologic image of the reservoir bttil tprovided by the 

geologist from all the information ho has gathered on the reservoir fracturation (this 
] 5 validation being performed by comparison with the well test data), 

reliably tes ttestinp the sensitivity of the hydraulic b ohaviou rb ehavior of the 

medium to the uncertainties on the geologic image of the fractured medium. 

A well-known modelling method consists iin s finely meshing the fracture network 
and the matrix while making no approximation concerning fluid exchanges between the 
20 two media. It is however difficult to implement because the often complex geometry of 
the spaces between the fractures makes it difficult to mcah th e m ond the meshing of the 
spaces, if^ any case, the number of meshes to be processed is finally often 
tromondoug very large . The complexity increases still further with a 3D mesh pattern. 
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Techniques for modelling fractured porous media are described in French patents 
2,757,947 and 2,757957 filed by the ftpplioon t assignee . The fefmeF technique of the first 
patent relates to determination of the equivalent fracture permeability of a fracture 
network in an underground multilayer medium from a known representation of thkthe 

5 network, allowing to aystomatioQlly oonnoo t svstemattc connection of fractured reservoir 
characterization models to double-porosity simulators in order to obtain more realistic 
modelling of a fractured underground geologic structure. The seeend technique of the 
second patent relates to simplified modelling of a porous heterogeneous geologic 
medium (such as a reservoir crossed through by an irregular network of fractures for 

10 example) in the form of a transposed or equivalent medium so that the transposed 
medium is equivalent to the original medium, according to a determined type of 
physical transfer function (known for the transposed medium). 

SUMMARY OF THE INVENTION 

The method according to the invention allows to simulat e simulation single-phase 
15 flows In a ftactured medium. 

The method can be implemented for example in the field of oil production by 
reservoir engineers in order to obtain reliable flow predictions. 

A well test is to be simulated in a permeable porous medium freseivoir'^ crossed 
through by a network of fractures much more conducting than the porous matrix. 
20 Fractured reservoirs constitute an extreme type of heterogeneous reservoirB comprising 
two very different media, a m atrix medium containing the most part of the oil in place 
and having a low permeability, and a fractured medium representins less than I % of 
the oil in place and highly conducting. The fractured medium ^tself cap be complex. 
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with different scries of fractures characterized bv their respective density, length, 
orientation, inclination and opening. 

During well testing, the flow rate conditions imnosed on the well lead the oil 
contained in the reservoir to flow towards the well. It is a sinple-phasc (the only mobile 
5 phase is the oil) and compressible flow fthe rock of the reservoir and the oil fluid arc 
commessibley 

For any elementary volume of the reservoir, the pressure of the oil contained in this 
volume is governed bv the following equation, if gravity is not taken into account : 

^,C^,^ = dH-.VP)-^Q (1) 
dt ft 

where : 

10 (b : represents the nore volume. 

Ct. the total comnressibilitv ffluid + rockV 

K, the permeabilitv of the rock, 
the viscosity of the fluid, 

O. the incoming flow, which is zero everywhere except in places where the well 
15 communicates with the reservoir, and 

P. the unknown pressure. 

In order to simulate a well test whatever the medium, this equation has to be solved 
in space and in time. Discretization of the reservoir (mesh pattern) is therefore 
performed and solution of the problem consists in finding the pressures of the meshes 
20 with time, itself discretized in a certain number of time intervals. 
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In the case of a fractured medium, if all the complexity of the fractuie network is to 
be taken into account it is necessary to have a mesh pattern that represents this network 
accurately. Furthermore, in such a medium, the permeabilities K are generally much 
higher in the fractures than in the matrix, and the flows are consequently fast in the 
5 fractures and slow in the matrix. 

The objoGt of tho method according to the invention i9 - te - <nodol models fluid flows 
in a fractured multilayer porous medium by taking aooount o fa ccounting for the real 
geometry of the fracture network and the local exchanges between the porous matrix 
and the fractures at each node of the network, and thus to allow oimulation o fs imulates 

10 the interactions between the pressure and flow rate variations in a well running across 
saidthe medium. It i s charactoriaod in that tho The fractured medium is discretized by 
means of a mesh pattern wherein the fractuie meshes are ©eHteedcentered on nodes at 
the various intersections of the fractures, each node being associated with a matrix 
volume, and the flows between each fracture mesh and the associated matrix volume are 

1 5 determined in a pseudosteady state. 

The matrix volume associated with each fracture mesh in each layer is for example 
d e limit e d defined by all of the points that are closer to the corresponding node than to 
noicJibourin g neighboring nodes. 

In order to determine the matrix volume associated with each fracture mesh, each 
20 fractured layer is for example discretized in pixels and the distance from each pixel to 
the closest fracture mesh is determined. 
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The value of the transmissivity for each ftacture mesh - matrix block pair is 
determined for example by considering that (he pressure varies linearly as a function of 
the distance from the point considered to the fracture mesh associated with the block. 

The method according to the invention greatly simplifies calculation of the 
5 exchanges between the matrix and the fractures. In fact^ there is no specific mesh 
pattern for the matrix, but each node of the fracture network is associated with a matrix 
volume, and the flow between each fracture node and its associated block is calculated 
in a pseudosteady state. The innovatio n inyention concerns calculation of the volume of 
each block and of proportionality coefficient (called transmissivity) that connects 
1 0 the matrix-fracture flow rate to the matrix-fracture pressure difference. 

The main advantage of the method according to the invention in relation to methods 
applied to a finely meshed real fracture network lies, as detailed in the description 
hereafter, in *ea considerable reduction in *ea number of meshes eegfee ^entered on 
fracture nodes used to solve equation (1). It can bo Qhoolcod that thoT he simulation rate 
15 saving in relation to prior fine-mesh methods is can be confirmed to be greater than the 
simple arithmetical ratio of the number of meshes used in both cases. 

BRIEF DESCRIPTION OF THE DRAWINGS 

Other features and advantages of the method according to the invention will be 
clear from reading the description hereafter of a non limitative realisation example, with 
20 reference to the accompanying drawings wherein : 

- Figuf&Fig. 1 shows an example of a 2D firacture mesh, P M Fracture Mesh . 

- FtgUF&Fig. 2 shows an example of a 2D matrix block, MB Matrix Block. 
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Figur e Fig. 3 shows an imposed flow rate variation curve F(t) in a well test, 

- Figia eFig. 4 shows an example of a fracture network for which the method 
according to the mvention leads to a fadtea teubstantial reduction in the number of 
meshes to be processed , and 

5 - Figs. 5 and 6 illustrate flow charts of the method of the present invention . 

DETAILED DESCRIPTf ON OF THE METHOi tfNVENTION 
1) Fracture network discretization 

The method comprises meshing the multilayer fractured reservoir modelled for 
example by means of the fractured reservoir modelling technique described notably in 
10 the aforementioned French p atent 2,75 7^,947, assuming that the fractures are 
substantially perpendicular to the layer boundaries. 

In Older to form this mesh pattern, all the fracture intersections or computing nodes 
GN are located (Fig.l) and fracture meshes are formed by associating a single matrix 
block MB (Fig.2^ with each node. For the purpose of 3D modelling, additional nodes 

15 have to be added in the n e ighbourin gn eighboring layers in order to account effor 
the fractures that run across several of th s m layers . The c ompurin g computed nodes thus 
defined constitute the eenfeecenter of the fracture meshes PM. The meshes are given 
boundaries such as the ends of the fractures on the one hand and the midpoint of the 
segments connecting the computing nodes on the other hand. Considering these 

20 boundaries imposed on the meshes, thetfthe volume ^ of the meshes can be calculated 
(seeSee relation 1). Calculation of the connections between fracture meshes 
(transmissivities) used for calculation of the inter-mesh flows can be carried out 
according to the method described in the aforementioned French patent Kft-2,757,947. 
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2) Matrix inediiim discretization 

Discretization of the matrix medium conoipto in ossigniH ^ssiens a rock volume to 
each fracture mesh. During dynamic simulation, exchanges between the matrix and the 
fractures wiH-beare calculated in the pseudosteady state (exchange flow proportional to 
5 the pressure difference) between each fracture mesh and Itsthe associated single matrix 
block. 

For calculation oftiwsettie matrix volumes, the problem is dealt with layer by layer, 
which amounts to disregarding the vertical flows in relation to the horizontal flows in 
the matrix. 

1 0 For a given layer, the blocks are defined as follows : 

- the block heights are equal and fixed to the height of the layer, 

- in the plane of the layer, the surface area of the block associated with a fracture 
mesh is all the points that are closer to ^tisthe fracture mesh than to another one. 

Physically, this definition implies that the boundaries at zero flow in the matrix are 
13 the equidistance lines between the fracture meshes. This approximation is valid if the 
n e ighbourin gn eighboring fractures are at very close pressures, which is true in the 
presence of highly conducting fractures in relation to the matrix. 

a) Block geometry 

In order to determine the geometry of the blocks in a given layer, a two-dimensional 
20 problem oonaisting in findin gfa as to solved which finds , for each fracture mesh^ the 
points that are closer to this mesh than to another one has to b e solv e d . This problem is 
solved by using for example, but preferably, the geometric method described in the 
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Other aforementioned French p atent 2,757,957, which allows, by discretizing the 
^ctuied layer into a set of pixels and by applying an image processing algorithm, to 
determine the distance from each pixel to the closest fractme. During the initiahzation 
phase of this algorithm, value 0 (zero distance) is assigned to the pixels belonging to a 
5 fracture and a high value is assigned to the others. Furthermore, if the number of the 
fracture mesh to which each fracture pixel belongs is given in this phase, the same 
algorithm finally allows to determine, for each pixel : 

- the distance from this pixel to the closest fracture mesh, 

- the number of the closest fracture mesh. 

10 All of the pixels having the same fracmre mesh number constitute the surface area 
of the matrix block associated with this mesh. Multiplying this surface area by the 
height of the layer allows to obtai n obtaining the volume of the block associated with the 
mesh. gjffUf eFig. 2 shows an example of a thus obtained 2D matrix block. 

By constmction, the sum of the surface areas of the blocks is equal to the total 
1 5 surface area of the layer, which guarantees that no volume is added to or taken from the 
layer. 

For each matrix block, the image processitig algorithm also gives the distance from 
each pixel of the block to the associated fracture mesh. This data is used to calculate the 
transmissivity between the fracture mesh and the matrix block. 

20 
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b) Matrix-fracture transmissivity 

On the assumption of a pseudosteady state where the exchange flow is considered 
to be proportional to the difference in the pressures of each one of them, wo hoi i 'O the 
following relation exists : 

(2) 

'ft 

5 where : 

Fmf is the matrix-fracture flow^ 

Tmf, the matrix-fracture transmissivity, 

|i, the viscosity, 

Pffl, the pressure of the matrix block, and 
1 0 Pr, the pressure of the associated fracture mesh. 

On the assumption of a pseudosteady state, the value of transmissivity is 
constant during simulation. This value is here calculated for each fracture mesh - matrix 
block pail. 

For this calculation, it is assumed that, in the matrix blocks the pressure varies 
15 linearly as a function of the distance from the point considered to the ftacture mesh 
associated with the block. The transmissivity is defined as follows : 

2\.HX 

where : 
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12 

1 is the length of the fractures in the fracture mesh, 
H, the height of the layer (and of the matrix block), 

the permeability of the matrix, and 
D, the distance from the fractures for which the pressure is the mean pressure of the 
5 block. 

1, H and K are known (factor 2 comes from the feet that the fracture has two faces). 
Calculation of D is made from the information provided by the image processing 
algorithm concerning the distances from the pixels to the closest fractures. In feet, 
according to the hypothesis of linear variation of the pressure as a function of the 
10 distance to the fracture, w e havo tfae relationship is : 

D = l.ldis).ds _(4) 
s 

where S is, in 2D, the surface area of the matrix block and d(s) is the function giving the 
distances between the points of the matrix block and the associated fracture mesh. 

In tenns of pixels, w e thu s s imply hav e t he relationship is : 

15 where N is the number of pixels of the matrix block and d„ the distance from pixel n to 
the fracture mesh. 

3) Well representation 

A well is represented by its geometry on the one hand and by the flow rate 
constraints imposed thereon with time on the other hand. 

20 
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a) Geometry 

A well oonaioto ofi s a series of connected segments that intersect the network 
fractures. The geometric representation of a well is therefore a 3D broken line. Some 
segments of this broken line communicate with the reservoir. The points of 
5 communication between the well and the network are the points of intersection of these 
segments communicating with the fractures. 

b) Flow rate constraints 

The flow rate constraints are imposed at the point where the well enters the 
fractured medium. ^^Phev The flow rates are entered by the user in the form of a step 
10 function giving the flow rate as a function of time. The flow rate change times of the 
well indicate the various periods to be simulated. 

For a conventional well test (Fig.3), a flow rate F(t) is imposed for a first period 
after which the well is closed (zero flow rate) and the pressure buildup is observed in 
the well. 

\S 4) Dynamic simulation 

During dynamic simulation, the simulated period of time is split up into time 
intervals whose duration ranges, for a well test, between 1 s (just after opening or 
closing of the well) and aem ein a range of hours. 

Passage from the time n (where all the pressure values are known) to the time n+1 
20 is achieved by solving, in all the meshes and blocks of the reservoir, the following 
equation which corresponds to relation 1 once discretized : 
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where : 

i, the number of the mesh or of the block, 
the number of the mesh or of the block next to i, 
5 P", the pressure at the time n, 
t", the time at the time n, 
Qi, the flow entering i, and 
Tij, the connection transmissivity between i and j. 

Apart from the pressures, the terms of this equation are known. The pore volumes 
10 of the fracture meshes and of the matrix blocks {^{) are known by means of the mesh 
pattern, the fracture-fracture and matrix-fracture transmissivities (Ty) are calculated as 
described above and flow rates (Qi) are zero everywhere except at the well-reservoir 
connections where they are imposed. 

a) ^liminatioii of the unknowns relative to the matrix 

15 Equation ^(7^ can be solved by inversion of a linear system where the unknowns 

are the pressures of the fracture meshes and the pressures of the matrix blocks at the 
timen+1. 

However, since each matrix mesh is only coimected to one matrix block and since 
the blocks are not connected to each other, it is possible to eliminate the unknowns of 
20 the blocks of the system whose size is then divided by 2. 

In fact, for a matrix block, equation is expressed as follows : 
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subscripts m and f denoting the matrix and the fracture respectively, and d™ being the 
total compressibility of the matrix. 

W&^ ^The i^lationship is Ciw - c,, + Sw.c^ + (l-Sw).Co, where Sw is the residual 
water saturation in the matrix, Cn, the compressibility of the matrix, Cw the 
5 compressibility of the water and Co the compressibility of the oil, and is the pore 
volume of the matrix block, irftr that is the volume of the block multiplied by the porosity 
of the matrix. 

We The relationship is deduce therefrom that : 

p^^i ^ _K — H- ^" - i^r' (8) 
with : 

10 By transferring this expression into the equation relative to the fracture mesh, we 
ebtei fithe relationshin is obtained ; 



(9) 



where i is the number of the fracture mesh considered, j the numbers of the fracture 
meshes connected to i, and m the number of the block associated with fracture mesh i. 

15 Furthermore, wo havo t he relationshin is : 
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where Ott is the total fracture compressibility : Crt = Cr + Co- 

Reduction of the number of unknowns by a factor 2 is a further advantage of the 
present method. The solution time for a linear system develops approximately like the 
square of the number of unknowns, and consequently elimination of the matrix 
5 unknowns leads to a considerable saving in the rate of simulation of a well test on a 
fractured network. 

h) Time interval solution algorithm 

The general algorithm for solving a time interval oonGioto iu foimhi gf orms a linear 
system corresponding to equation (9) relative to any fracture mesh i, in solving this 
1 0 linear system and in updating the unknowns. 

The stages of this algorithm are as follows : 

- Forming the linear system with calculation of the terms of equation (9) 
Calculation of the fracture medium contribution 
accumulation term (A;) 
1 5 connection between fracture meshes (Ti/n) 

well boundary conditions (Qj) 

Calculation of the matrix medium contribution 
accumulation term (A„) 
matrix-fracture connection term (Um) 

20 - Solving the linear system (so as to find the fracture pressure values at the time n+l) 
by means of the iterative conjugate gradient algorithm (COS) 
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- Updating with calculation of the matrix block pressures using equation (8). and time 
interval management. 

c) Time interval management 

During well test simulation, the time intervals are generally very short after opening 
5 or closing of the well, because the pressure variations are high at these times. Later, the 
time intervals can be longer when the pressure variations are lower. 

Management of the time intervals oonoiDto in oonnootin g connects the duration of the 
time interval n+1 to the maximum pressure variation observed during time interval n. 
The user roust therefore provide the following data : 

1 0 Initial time interval : dtb,i 

Minimum time interval ; dtmin 

Maximum time interval : dt„j„ 

Lower limit of pressure variation : dpn,i„ 

Upper limit of pressure variation ; dprau 
15 Time interval increase fector : rdt + (>1) 

Time interval reduction factor : rdt - (<1). 

From these values, management of the time interval is performed as follows : 

At the time t,,, the H iatoion t maioritv of the pressure variations in the fracture meshes 
and the matrix blocks during time interval n is calculated (i.e. between times f ' and f), 
20 According to the value of this pressure variation dp, the time interval is either : 

- increased by a factor rdt+ if dp < dp^ri,,, 

- decreased by a &ctor rdt- if dp > dpnw«* 
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unchanged in the other cases. 

After each well flow rate condition change, the time interval is re-initialized to the 
value dtfa^. Furthermore, an optimization is performed at the end of the simulation 
period ; if tf is the time to be simulated to reach the end of the period and dt the planned 
5 time interval, the time interval selected is min(tffi,dt). 

5) Simulation input data 

a) Fracture network 

The geometry of the jfiacture network and the attributes of the fractures 
(conductivity, opening) are given in the form of a file as in the method described in the 
1 0 aforementioned French patent ¥Rr-2 ,75 7,947. 

b) Petro physics 

The petrophysical data to be provided are as follows (the unit is given between 
parentheses) : 

- compressibility of the fracture network (Cf in bar') 
1 5 - compressibility of the matrix (c„ in bar') 

■ matrix permeability (K in mD) 

- matrix porosity (4 . dimcnsionlo ss ^ for computing the pore volume . 

The data in question are considered to be homogeneous on the ftactured medium 
considered. 

20 c) Fluids 
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The data relative tg the fluids contained in the fractured medium concern the oil 
(mobile) and the water (immobile) that is always present in residual amounts in the 
matrix : 

- residual water saturation in the matrix (Sw dimensionless) 
5 - compressibility of the water (c« in bar ') 

- compressibility of the oil (c„ in bar') 
viscosity of the oil in cpo). 

d) Weil 

The data relative to the well are i 
10 - its geometry, in the form of a series of connected segments 

- the imposed flow rates, in the form of a curve giving the imposed flow rate as a 
function of time. 
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e) Time management 

The values to be provided are the values already defined in the chapter relative to 
the time interval management : dti,i, dUu, dtw,, dp^^, dp^, rdt+ and rdt-. 

Comparative meshing example 

5 The following comparative example illustrates the meshing simplification and the 

correlative saving in flow calculating time that can be obtained by implementing the 
method. By way of comparison, it has been observed that the fractm^d zone shown in 
Fig.4, that normally would have been meshed with 80,000 meshes with conventional 
meshing techniques, can be meshed here with about 4000 meshes only. 

10 Fig. 5 illustrates a flow chart of the overall processing of the invention. Processing 
starts at step 10. where the fracture network is defined as discussed above under 
fracture network definition. Processing proceeds to step 12 where the well is defined as 
discussed above u nder well representation. Processing proceeds to step 14 where the 
matrix medium is defined as discussed above with respect to matrix medium definition. 

15 Processing proceeds t o step 1 6 where computation of the transmissivities occurs as 
discussed aboveunder matrbc-fracturc transmissivitv. elsewhere with respect to the 
transmissitivity Tij and as further discussed with respect to a linear system as set forth in 
equation 9. Processing proceeds to step 18 where initialization occurs as discussed 
above with reference to fra cture networV and matrix medium discretization. Finallv. 

20 processing proceeds to step 20 where a dynamic simulation loop occurs involvmg time 
deffnitign. 

Fig. 6 illustrates a flow chart of the dynamic simulation loon 20. Processing in the 
dynamic simulation loop 20 begins with step 100 where building of the linear system 
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using data at time steps n. occurs as discussed above regarding time interval solution 
algorithm. Processing proceeds to step 102 where resolution of the linear system is step 
n+1 occurs as discussed under time interval solution algorithm, firnnessing proceeds to 
step 104 where time step management and update of well rate occurs as discussed under 

5 time interval management. Processing proceeds to decision point 106 where 

determination is made if an end of simulation has occurred. If an end of simulation is 
not sensed at decision point 106 as described abovejinder time interval management . 
the processing returns to point 100. If the end of simulation is sensed at decision point 
106. an output of the simulation which represents the pressure evolution of the well are 

10 produced as indicated at step 108. 
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